home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Languages / RLaB 1.18c / testmatrix / chebvand.r < prev    next >
Text File  |  1994-12-20  |  2KB  |  65 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Vandermonde-like matrix for the Chebyshev polynomials.
  4.  
  5. // Syntax:    C = chebvand ( P )
  6. //        C = chebvand ( M , P )
  7.  
  8. // Description:
  9.  
  10. //    C = chebvand(P), where P is a vector, produces the (primal)
  11. //    Chebyshev Vandermonde matrix based on the points P, i.e.,
  12. //    C(i,j) = T_{i-1}(P(j)), where T_{i-1} is the Chebyshev
  13. //    polynomial of degree i-1.
  14.  
  15. //      chebvand(M,P) is a rectangular version of chebvand(P) with M
  16. //    rows. Special case: If P is a scalar then P equally spaced
  17. //    points on [0,1] are used.
  18.  
  19. //      Reference:
  20. //        N.J. Higham, Stability analysis of algorithms for solving confluent
  21. //        Vandermonde-like systems, SIAM J. Matrix Anal. Appl., 11 (1990),
  22. //        pp. 23-41.
  23.  
  24. //    This file is a translation of chebvand.m from version 2.0 of
  25. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  26. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  27.  
  28. // Dependencies
  29.    rfile seqa
  30.  
  31. //-------------------------------------------------------------------//
  32.  
  33. chebvand = function ( m , p )
  34. {
  35.   local (m, p)
  36.  
  37.   if (!exist (p)) 
  38.   { 
  39.     p = m;
  40.   }
  41.  
  42.   n = max(size(p));
  43.  
  44.   //  Handle scalar p.
  45.   if (n == 1)
  46.   {
  47.     n = p;
  48.     p = seqa(0,1,n);
  49.   }
  50.  
  51.   if (!exist (_p)) { m = n; }
  52.  
  53.   p = p[:].';        // Ensure p is a row vector.
  54.   C = ones(m,n);
  55.   if (m == 1) { return C; }
  56.   C[2;] = p;
  57.   //      Use Chebyshev polynomial recurrence.
  58.   for (i in 3:m)
  59.   {
  60.     C[i;] = 2.*p.*C[i-1;] - C[i-2;];
  61.   }
  62.  
  63.   return C;
  64. };
  65.